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O ' Abstract 

c/j - 

We perform an extrapolative analysis of "fast-growth" free-energy-difference (AF) 

estimates of a computer- modeled, fully-solvated ethane^methanol transformation. 

The results suggest that extrapolation can greatly reduce the systematic error in AF 

estimated from a small number of very fast switches. Our extrapolation procedure uses 

block-averages of finite-data estimates, and appears to be particularly useful for broad, 

non-Gaussian distributions of data which produce substantial systematic errors with 

insufficient data. In every tested case, the extrapolative results were better than direct 

^ , estimates. 
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1 Introduction 

Relative free energy computations have long been of interest, and biological applications 
promise to be of particular importance |j], 0, j!|. As examples, it would be desirable to 
accurately and rapidly estimate free energy changes resulting from the opening of an ion 
channel, the binding of a ligand, and alchemical mutation among a series of protein ligands. 
Ligands might include potential drug compounds or varying sequences of nucleic acids (RNA 
and DNA). Strategies for computing free energy differences date back to Kirkwood M and 
Zwanzig || who pioneered thermodynamic integration and free-energy perturbation strate- 
gies. Many computational strategies have since been developed for molecular systems (e.g., 

BBBD- 

"Fast-growth" methods || 0, H, [§ HU HH' 0' HH are ^ ne f° cus °f the present paper. 



The impetus for these approaches comes from the work of Reinhardt, Hunter, and coworkers 
H [7| who recognized that computations could readily employ a microscopic analog of the 
inequality between work and free energy. The principle is readily illustrated in an "alchem- 
ical" context where one wishes to compute the free energy difference between two systems 
described by different potential energy functions, Uq and Ui, and parameterized by the 
switching variable A according to an extended potential function: 

C/(x; A) = f/ (x) + A[ ^(x) - f/ (x) ] , < A < 1 , (1) 

where x is a set of configurational coordinates. If one performs a series of rapid "switches" 
(described below) between the two systems using an amount of work W in each switch, the 
free energy difference is bounded according to [1_L Q 

- (Wi_+o) < AF Q ^ < (Wo-n) , (2) 

where the (• • •) brackets indicate averages over many switches starting from equilibrium 
ensembles of either start (A = 0) or end (A = 1) systems. (The distinct "systems" could 
also describe different conformations of a single system constrained to distinct values of a 
reaction coordinate.) 

The potentially rapid, non-equilibrium events used to compute (W) in Eq. (g) thus 
provide a computational estimate of the equilibrium quantity AF. However, the bounds will 
not be tight unless the switches are sufficiently slow, offsetting some of the computational 
savings. 

Subsequent work by Jarzynski ||, [K| sidesteps, at least in principle, some of the limita- 
tions by permitting direct computation of AF from a single set of rapid switches, via the 
simple, exact relation, 

e -AF/k B T = ( e -W/k B T^ (3) 

However, estimates for AF generated using Eq. (0) are highly sensitive to small values of W 
and significant errors can arise when the width of the distribution of W values exceeds k^T 



TT| , P~3|| . Hummer's recent work with a small molecular system concluded that little, if any, 
advantage was gained from the fast-switching approach [[13|] . 

In the past improvements have been sought in the procedure for generating a set of work 
values {Wi, W2, . . .} to be analyzed according to Eq. (H) or @. In particular, one can switch 
between A = and 1 along arbitrary paths, perhaps using more than one switching pa- 
rameter as initially discussed by Reinhardt and coworkers for the fast-switching approach by 



3, H]. Subsequent exploration of optimal switching paths has been pursued by many workers 
15| , D, [12], 0. In fact, the exploration of different paths in alchemical free-energy compu- 



tations pre-dates the fast-switching approach, and was pursued in free-energy-perturbation 
and thermodynamic integration efforts 



e.g. 



16, 171 18 



The present study, by contrast, attempts to optimize the use of the data {W\, W2, ■ ■ ■} 
which has already been generated, by using a combination of block-averaging and extrapo- 
lation. This additional statistical analysis is needed to bypass the systematic error inherent 
in finite data samples 0, ^, |Hj]. Fig. p] illustrates the basic points. The running averages 
(solid lines) based on Eq. (Q) exhibit erratic behavior, and it is essentially impossible to judge 
from these whether the computation has converged to an answer. However, the same data 
considered in block-averages (error bars) is well-behaved and, as seen below, well-defined. 
Only the block-averages could be considered for extrapolation to the "infinite-data" limit. 
Cases of insufficient data requiring extrapolation are of great interest because the size of 
biomolecular systems often makes relative free energy estimates extremely costly. 
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Figure 1: Running and block averages for forward and reverse switching. The evolving 
estimates for the free energy difference are plotted vs. the number of switches, N = N switch . 
The running averages [based on Eq. ([!]); solid lines] exhibit non-monotonic, rise-and-drop 
behavior, while the block averages (defined in Sec. 0; error bars) are monotonic and smooth. 
Each block-average data point was computed using all 10,000 work (W) values. The error 
bars are twice the standard error of the mean (see Sec. |3|), and represent roughly 90% 
confidence intervals [EI] . Data from switches of 20 A steps. 



Following Jorgensen and Ravimohan [22| and Jarque and Tidor |]], we examine alchem- 
ical mutations between methanol and ethane in explicit water solvent. The authors are 
unaware of any previous application of Jarzynski's relation to alchemical transformations in 
a molecular system, although Hummer performed a methodical study of the inter-methane 
distance dependence of the free energy |13| . 

Our results indicate that the combined use of block-averaging and extrapolation is very 
promising and warrants additional investigation. The approach produces successful and 



reasonably reliable relative-free-energy estimates even from very fast switches of only one 
or two steps, which generate extremely broad, highly-non-Gaussian distributions of work 
values. In every case we examined, extrapolation of the data yielded a better estimate than 
direct averaging alone. 

In outline, this Letter is organized as follows: Sec. |2| briefly describes "fast-growth" 
computations and gives simulation details. In Sec. |3] we define the block averages, and the 
extrapolation procedure is discussed in Sec. [|. In Sec. |5]we summarize our results and discuss 
future work, including potential applications of the approach to large biomolecular systems. 
We also discuss implications for other approaches to free energy calculations. 



2 Alchemical Free Energy Calculations 

This section fills in some details regarding the theory governing an "alchemical" free en- 
ergy change and its implementation using a rapid-switching strategy. Alchemical changes 
are transformations between Hamiltonians which describe different molecules; molecular iso- 
merization is mathematically analogous but not considered here. The free energy difference 
between the two states is formally given by the ratio of the partition functions according to 

exp (-AiWfcT) = J fd ^ e _ Uo{x)/kBT ■ (4) 

Jarzynski's relation (|J) is derived from this definition. 

Free Energy Perturbation (FEP) 

The so-called free-energy-perturbation (FEP) procedure for computing relative free energies 
||, ^ [23|] is a well-established method for molecular systems |T7|, [18|, [I], 0, |3| which we use 
as a benchmark for understanding systematic errors. FEP computations entail a number of 
equilibrium simulations performed at a set of fixed values of A; for example, our FEP result 
quoted in Sec. [| uses simulations at A = 0.1,0.2, . . . ,0.9. The total free energy change is 
estimated as the sum of the incremental changes, which are computed based on the analog 
of Eq. (Q) involving (exp (— AW/k^T)), where AW is the work or energy difference between 
configurations at different A values. 

Fast-Growth Procedure 

Fast-growth algorithms have been discussed in detail elsewhere (e.g., [^, |TU|, [TTj, [T3J), so we 
merely sketch the approach. The general procedure for computing a "fast-growth" free energy 
difference - - via Eq. (0) rather than (f|) - - begins with the generation of an equilibrium 
ensemble of starting (A = 0) configurations, perhaps by molecular dynamics simulation as 
is done here. One proceeds by (i) choosing a configuration from the equilibrium ensemble, 
(ii) incrementing the potential energy function (|l|) to a new, greater value of A (keeping the 
configuration fixed) and (iii) relaxing the system at the new A value. Steps (ii) and (iii) 
are repeated until a value A < 1 is reached. In our implementation, the A increments in 
(ii) are uniform and the relaxation stage (iii) consists of a single molecular dynamics (MD) 



"relaxation" step, following the "fast-growth" convention [|5], 1I|. A uniform increment of 
AA = 0.05, for instance, corresponds to 20 "A steps" and would require 19 MD steps, as 
none is necessary at A = 1. 

The work for any such switch is computed based only on the potential energy increments 
and not the relaxation dynamics. Thus, if xf n denotes the final configuration of the system 
after it is relaxed at the ith value A;, the work calculated from 



»u 



w = E M*Ei> A *) - u (*£i> M] > ( 5 ) 

where the same configuration is evaluated at two different A values. Finally, to evaluate the 
averages in Eqs. (0) and ([3]), one uses additional members of the A = equilibrium ensemble 
to generate subsequent values of W — starting from step (i), above. 



Met hanol^ Ethane Model and Simulation 

Simulations of the methanol«-»ethane "transmutation" were performed within the CHARMM 
molecular dynamics package. Both methanol and ethane were modeled in the united-atom 
picture: methanol was represented as a three-atom (C,0,H) molecule and ethane as a two- 
atom (C,C) molecule. The solvent used 125 TIP3 water molecules (for both A = and 1) in a 
periodically replicated box of (15.6 A) 3 . To facilitate comparison with earlier studies, electro- 
statics and van der Waals interactions were both shifted to zero at a cutoff of 8 A. Molecular 
dynamics steps (performed at fixed A values) used the leapfrog Verlet algorithm. The same 
simulation procedure and parameters were used for free energy perturbation calculations. 



3 Block Averaging 



While block-averaging is straightforward, its repeated application for growing block sizes 
to a non-linear transformation - - such as taking the log of an average of exponentials in 
Jarzynski's relation (|3|) - - turns out to yield rich, well-behaved data: see Fig. [[]. The 
procedure and some implications are discussed now. 

We construct block averages [E_0, E_4[ |K], |25| from a set of, say, N tot work values 
{V^i, W2, • • • , Wiv tot } by applying Jarzynski's relation ([5]) to a series of blocks, each containing 
iV = iVgwitch values. More specifically, we define the iV-block- averaged estimate for the free 
energy as 

iWJV 
AFn = W~ E -^ T1 °S (e- W/kBT )N, n = (f N ) , (6) 



n=l 



where the individual block averages are defined by 
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= (n-l)7V+l 



-Wi/k B T _ f 
e — j7V,n 



(7) 



The ratio N tot /N denotes the largest integer less than or equal to the literal fraction, and 
is never less than 30 in our analysis. Because of potential correlations in the sequence 
{Wi, W2, ■ ■ ■ , Wjvtot} we randomly re-sort the values prior to computing the block results 
presented here. We note that larger block sizes, N, could be considered with a bootstrap 
P^j or subsampling |^5] analysis. 

The true free energy difference of Eq. (^) is AF = Ai 7 ^, and the other limit gives 
the average work, (W) = Ai*\: see Fig. |I]. In general, a finite value of N indicates that 
the average in Eq. (|7|) is performed from a poor sample of the W distribution, with N 
determining how much of the tails of the distribution are included in the average. However, 
the averaging of these poor samples in Eq. (||) yields a well-defined descriptor of the finite- N 
statistics. In the present case, the Boltzmann-factor form ensures monotonic behavior, so 
that 

AF N+1 < AF N , (8) 



the essence of which was noted by Jarzynski JTOf ; see also |2(J . The usual relation between 
the average work and free energy (El) is simply a weaker case of the more general inequality 



The uncertainty in the finite- N free energy values, 5AF, is estimated by twice the stan- 
dard error of the mean, 



(SAF* ^ 



Ntot/N 



'N 



(iVtot/iV) 2 



J2 (fN,n-(f N )f 



n=l 



(9) 



which gives roughly a 90% confidence interval [21]. This is the quantity used to compute 
error bars and uncertainties. 



4 Extrapolation 

While extrapolation and data-fitting are something of black arts, one can hope to derive 
meaningful information with a careful error analysis [p6fl . Here we discuss some simple, 
intuitively appealing schemes for extrapolating the block-averaged, finite-data free energies 
(H) to the limit of infinite data. The motivation for our approach is the analysis of finite-sue 



effects in spin systems [E7], [28 



Inspection of the data on a linear scale, such as Fig. [I], and in logarithmic plots suggests 
the simplest fit might be to a power law, 

AF N = AF OQ + a 1 (l/N) a \ (10) 

A natural, related form considers a power series 

ktnax 

AF N = AF 00 +J2h(l/N) k ^, (11) 

fe=i 

where the parameter (3\ can be chosen from a fit or some other way, such as by examining the 
leading 1/N behavior. Our work with the form (|TT| ) uses three parameters with k max = 2, 
except where noted, and the fixed exponent /3i = 0.266 chosen empirically, but based on 
some of the values fitted for a± in Eq. (|Il]). Naturally other exponents and polynomial 
degrees could be used. 

One drawback to these forms is clear: if the data do not include the leading 1/N behavior 
and the "distance" to extrapolate is great (from 1/N = to the first data point; see Fig. |2|), 
the fits will not have good extrapolative power. We anticipate that an analytic understanding 
of the behavior of AF^(N) for model systems, to be pursued in future work, will shed light 
on extrapolation forms and methods. 



5 Results 

We now present estimates for the free energy difference of the methanol— methane trans- 
formation, based on the block-averaging and extrapolation presented in the previous two 
sections. Our focus is the methanol-to-ethane direction of the transformation because it is 
more challenging and so presumably a better model for larger systems. 

The basic results are surprising and exciting. First, successful extrapolation to reason- 
ably accurate free energy values does appear to be possible in the methanol— ►ethane system. 
Moreover, for fixed amounts of computer time, the extrapolated estimates appear to be con- 
siderably better than standard fast-growth values, and can avoid errors of several kcal/mole 
resulting from insufficient data. If borne out for other systems, the ability to make estimates 
from a relatively small number of very rapid switches would mean dramatic efficiency gains. 
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Figure 2: Extrapolation of the free energy estimate for a fully-solvated methanol^ethane trans- 
formation. Finite-switch averages A-F/y are plotted as a function of the number of switches per 
average iV = N sw - ltc ^ raised to a "scaling" power. The block averages, fit, and direct estimate were 
all computed from the same data, while the free-energy-perturbation (FEP) value was generated 
from an independent, substantially longer calculation. The symbols roughly indicate the sizes of 
the error bars. The block-averaging is described in Sec. |j[ the extrapolative fitting in Sec. ||, and 
the direct estimate averages the same data according to Eq. (JJ). The data are from 10 3 switches 
of 20 A steps each, and note that ksT ~ 0.6 kcal/mole. 



Fig. H shows a sample extrapolation, based on Eq. ([LTD , f° r 10 3 switches of 20 A steps 
each. Note that the un-extrapolated, "direct" free-energy estimate -- based on application 
of Eq. @ to the same data -- exceeds both the extrapolated value and the reliable FEP 
estimate by 7 kcal/mole~ llksT. Thus, with a limited amount of data, extrapolation of the 
block-averaged values yields a much better estimate. 



Table 1: Extrapolated estimates for the solvated methanol— methane free energy difference. 
The estimates AF cst and uncertainties are given in units of kcal/mole, and may be com- 
pared to the free-energy-perturbation estimate of 5.3 ±0.16 kcal/mole. Direct estimates 
are computed from Eq. (|3|) and power-series extrapolations from (|TTD , using identical data. 
The uncertainties in the extrapolations are discussed in Sec. [|. The bracketed values give 
the differences between the direct AF estimates and the more costly FEP estimate in the 
first row, and hence measure the accuracy of the former. The quantity "A steps" indicates 
the number of increments in the alchemical coordinate: see Sec. |l|. "Total Steps" gives the 
number of MD steps excluding those for generating the equilibrium ensemble at A = 0. 



Method 


A^pcst 


Uncert 'y 


A steps 


Tot. 


Steps 


Direct 


7.37 


[2.1] 


200 


2 


10 6 


Extrapolation 


5.68 


1.66 


200 


2 


10 6 


Direct 


11.2 


[5.9] 


200 


2 


10 5 


Direct 


8.50 


[3.2] 


20 


2 


10 5 


Extrapolation 


4.93 


0.960 


200 


2 


10 5 


Extrapolation 


6.41 


1.21 


20 


2 


10 5 


Direct 


12.7 


[7.4] 


20 


2 


10 4 


Direct 


8.68 


[3.4] 


2 


1 


10 4 


Extrapolation 


4.87 


1.08 


20 


2 


10 4 


Extrapolation 


7.85* 


1.46* 


2 


1 


10 4 



* These values change to 5.03 and 1.85 for k r , 
goodness-of-fit measure. 



3 in Eq. ([IT]) with a substantially improved 
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Table [I] presents quantitative results for the methanol-to-ethane transformation. The ex- 
trapolations are uniformly superior to the direct estimates for any fixed amount of computer 
time and consistently avoid errors on the order of several kcal/mole (where 1 kcal/mole 
~ 1.6 ksT) for smaller amounts of data. Total computer times for the tabulated results 
range from 2 nsec. (10 4 switches of 200 steps) down to just 10 psec (10 4 2-step switches) of 
non-equilibrium molecular dynamics simulation. The "Total Steps" column does not include 
the computer time expended on generating an equilibrium ensemble at A = because it is 
unlikely that one would investigate the transmutation of a system which has not already 
been subjected to an equilibrium study. 

We estimated upper and lower bounds simply by extrapolating, independently, from 
the sets of upper and lower limits of the confidence intervals; recall Eq. @. Statistical 
uncertainties were not given for the direct estimates because the systematic error is clearly 
more significant than the statistical: the bracketed deviations in Table p] indicate the direct 
estimates differ dramatically from the free-energy perturbation (FEP) result. Recall that 
the FEP approach was outlined in Sec. |[ 

The power of the extrapolative approach is underscored by the challenging character of 
the distributions of work values under consideration. The distributions are all quite broad 
and asymmetric: standard deviations range from 12 kcal/mole ~ 20 ksT (for the 200-step 
switches) to 24 kcal/mole ~ 38 ksT (2 steps), and third moments range from 72% of the 
standard deviation (200 steps) to 100% (2 steps). Thus, although all of the tabulated 
simulations involve very rapid switches - - of less than 1 psec. of molecular dynamics time 
per switch — the substantial differences in the distributions indicate that the data sets are 
quite distinct. We also noted a degree of robustness in trials with related but different 
forms and exponents (3\ (results not shown) which typically yielded consistent, if slightly 
inaccurate, results across data sets from widely disparate numbers of A steps — and hence 
disparate computer times and work distributions. 

Despite the success of the fitting form used here, superior extrapolations may be possible. 

The forms employed here (see Sec. |]) are empirical, so a theoretical basis should provide 

additional insight. Lacking that, a more systematic exploration of the implicit parameters 

- the exponent (3\ in (|Tl~D, the minimum number of switches per block, and the degree of 

the fitting polynomial — would also be valuable. 

Another interesting trend illustrated in the data of Table |I] is that for a fixed amount 
of computer time, direct estimates using fewer A steps appear to give better results. The 
statistical errors (data not shown) are also better for direct estimates using more rapid 
switches. 
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6 Summary and Discussion 

We have performed and analyzed extrapolative free energy estimates based on "fast-growth" 
alchemical simulations of a fully solvated methane^-^ethanol transformation. The results of 
Table [I] suggest that the combined use of block-averaging (Sec. Q) and extrapolation (Sec. |) 
permits accurate estimates from a relatively small amount of data which — when analyzed 
using the standard "direct" method - - leads to unacceptably large systematic errors of 
several kcal/mole. Extrapolated results, for our system, were always better than standard, 
direct estimates. The approach also appears to be fairly robust, in that good results are 
achieved over ranges both of overall computer time and of alchemical switching speeds. 

Our work builds on that of Wood et al, who perceptively proposed a first-order estimate 
of the systematic errors due to finite samples of data p0 |. The present method, however, is 



not limited to narrow work (energy-change) distributions as noted in Sec. [5]. 

This Letter describes an initial exploration of a potentially important approach, and 
a number of important issues and questions deserve further exploration. To name a few: 
(i) undoubtedly, simultaneous fits of forward (A = — > 1) and reverse switching data will 
provide more reliable free energy estimates; (ii) we have not performed a quantitative analysis 
of the efficiency, both by comparison to standard "fast-growth" approaches as well as to free- 
energy-perturbation estimates; (iii) how does the extrapolation approach generalize to larger 
biomolecular systems? (iv) how universal are the behaviors of the finite-data estimates, 
AF/v, considered in the extrapolation? (v) can theoretical scrutiny of simple models and 
distributions clarify the extrapolative procedure? The ideas discussed here may also apply, 
with suitable modifications, to perturbative calculations. 

We have discussed methods for analyzing data from fast- switching simulations, but have 
not broached the possibilities for improved sampling of data. There appear to be a number 
of promising, unexplored avenues. Instead of using a uniform alchemical increment AA, 
for example, one could adjust increments to ensure relatively constant work increments, 
following the example of perturbative calculations fl6], [17|, HH, |29| ; this approach could also 
be adapted for higher- dimensional alchemical coordinates already considered by others |], 0, 
|30| , §]. Improved sampling efficiency may also result from biasing the "relaxational," fixed-A 
dynamics to favor states with smaller work increments. 

Finally, we note that the relationship between the approach described here and estab- 
lished statistical methods needs to be elucidated. Elements of our approach, particularly 
the construction of "finite-data" block averages, clearly have been considered in "bootstrap" 
23J and "subsampling" EH] statistical approaches. Nevertheless, the authors are not aware 



of a similar practical -- if ad hoc -- technique for extrapolation to the infinite-data limit 
like that presented here. 
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